## load libraries for spatial analysis & plotting
library(sf)
library(spdep)
library(tidyverse)
library(classInt)
## read districts
# districts <- st_read("/home/kai/Work/data/gis/vg250_2014-01-01.utm32w.shape.ebenen-stand-januar-2014/vg250_ebenen/VG250_KRS.shp")
# This is copyrighted. You need to download the shapefile from here: 
# https://gdz.bkg.bund.de/index.php/default/verwaltungsgebiete-1-250-000-stand-01-01-vg250-01-01.html

# Two counties were merged
# First, change two old GKZ to one new
# change from factor to character
districts$AGS <- as.character(districts$AGS)
districts$GEN <- as.character(districts$GEN)
districts$AGS[districts$AGS == "03152"] <- "03159"
districts$AGS[districts$AGS == "03156"] <- "03159"
# Now, change the name
districts$GEN[districts$AGS=="03159"] <-"Göttingen (neu)"
# Do da merge
districts <- districts %>% group_by(AGS,GEN)  %>% summarise(geometry = sf::st_union(geometry)) %>% ungroup() 

# Flage up the new district to make sure that it has worked 
districts$neu <-0
districts$neu[districts$AGS=="03159"] <-1
ggplot(districts) + geom_sf(aes(fill=neu))

## read summaries from score
scoredata <- haven::read_dta("districts.dta")
## GKZ needs to be a string in both data frames
scoredata$GKZ <- scoredata$district
districts$GKZ <- as.character(districts$AGS)
## merge!
mergeddistricts <-inner_join(districts,scoredata,by="GKZ")

# Prepare for calculating Moran's I etc. 
# Find queen neighbours
neighbours <- poly2nb(mergeddistricts)
# Define weight matrix 
listw <-  nb2listw(neighbours)
# Calculate global Moran's Rs
populism.Moran <-  moran.test(mergeddistricts$populism, listw)

# save Moran to file 
sink("raw-populism-global-moran.txt")
cat(round(as.numeric(populism.Moran$estimate[1]),digits=3))
sink()


ggplot(mergeddistricts, aes(x=populism)) + geom_density() + xlab("Populism")
ggsave("populism-raw-distribution.pdf",units="mm",width=150,height=50)


## plot Islamophobia
## find quantiles
breaks_qt <- classIntervals(mergeddistricts$islamophobia, n = 7, style = "quantile")

mergeddistricts <- mutate(mergeddistricts, islamophobia_cat = cut(islamophobia, breaks_qt$brks,include.lowest=TRUE))

ggplot(mergeddistricts) + geom_sf(aes(fill=islamophobia_cat)) + scale_fill_brewer(palette = "YlOrRd",name="Islamophobia") + theme(plot.margin=grid::unit(c(0,0,0,0), "mm"))
ggsave("islamophobia.pdf",units="mm",width=150,height=150) 


islamophobia.Moran <-  moran.test(mergeddistricts$islamophobia, listw)

# save Moran to file 
sink("raw-islamophobia-global-moran.txt")
cat(round(as.numeric(islamophobia.Moran$estimate[1]),digits=3))
sink()


ggplot(mergeddistricts, aes(x=islamophobia)) + geom_density() + xlab("Islamophobia")
ggsave("islamophobia-raw-distribution.pdf",units="mm",width=150,height=50)



# plot populism 
breaks_qt <- classIntervals(mergeddistricts$populism, n = 7, style = "quantile")
mergeddistricts <- mutate(mergeddistricts, populism_cat = cut(populism, breaks_qt$brks,include.lowest=TRUE))
ggplot(mergeddistricts) + geom_sf(aes(fill=populism_cat)) + scale_fill_brewer(palette = "YlOrRd",name="Populism") + theme(plot.margin=grid::unit(c(0,0,0,0), "mm"))
ggsave("populism.pdf",units="mm",width=150,height=150)

# Plot cultural threat perceptions 
breaks_qt <- classIntervals(mergeddistricts$cultthreat, n = 7, style = "quantile")
mergeddistricts <- mutate(mergeddistricts, cultthreat_cat = cut(cultthreat, breaks_qt$brks,include.lowest=TRUE))
ggplot(mergeddistricts) + geom_sf(aes(fill=cultthreat_cat)) + scale_fill_brewer(palette = "YlOrRd",name="Cultural Threat") + theme(plot.margin=grid::unit(c(0,0,0,0), "mm"))
ggsave("cultthreat.pdf",units="mm",width=150,height=150)

cultthreat.Moran <-  moran.test(mergeddistricts$cultthreat, listw)

sink("raw-cultthreat-global-moran.txt")
cat(round(as.numeric(cultthreat.Moran$estimate[1]),digits=3))
sink()

ggplot(mergeddistricts, aes(x=cultthreat)) + geom_density() + xlab("Cultural Threat")
ggsave("cultthreat-raw-distribution.pdf",units="mm",width=150,height=50)


# Plot RWA
breaks_qt <- classIntervals(mergeddistricts$rwa, n = 7, style = "quantile")
mergeddistricts <- mutate(mergeddistricts, rwa_cat = cut(rwa, breaks_qt$brks,include.lowest=TRUE))
ggplot(mergeddistricts) + geom_sf(aes(fill=rwa_cat)) + scale_fill_brewer(palette = "YlOrRd",name="Authoritarianism") + theme(plot.margin=grid::unit(c(0,0,0,0), "mm"))
ggsave("rwa.pdf",units="mm",width=150,height=150)

rwa.Moran <-  moran.test(mergeddistricts$rwa, listw)

sink("raw-rwa-global-moran.txt")
cat(round(as.numeric(rwa.Moran$estimate[1]),digits=3))
sink()

ggplot(mergeddistricts, aes(x=rwa)) + geom_density() + xlab("Authoritarian submission/aggression")
ggsave("rwa-raw-distribution.pdf",units="mm",width=150,height=50)


# Perception of region as excluded 
breaks_qt <- classIntervals(mergeddistricts$regexclusion, n = 7, style = "quantile")
mergeddistricts <- mutate(mergeddistricts, regexclusion_cat = cut(regexclusion, breaks_qt$brks,include.lowest=TRUE))
ggplot(mergeddistricts) + geom_sf(aes(fill=regexclusion_cat)) + scale_fill_brewer(palette = "YlOrRd",name="Place resentment") + theme(plot.margin=grid::unit(c(0,0,0,0), "mm"))
ggsave("regexclusion.pdf",units="mm",width=150,height=150)

# We need moran for regexclusion, too
regexclusion.Moran <-  moran.test(mergeddistricts$regexclusion, listw)

sink("raw-regexclusion-global-moran.txt")
cat(round(as.numeric(regexclusion.Moran$estimate[1]),digits=3))
sink()

# And finally localism 
breaks_qt <- classIntervals(mergeddistricts$localism, n = 7, style = "quantile")
mergeddistricts <- mutate(mergeddistricts, localism_cat = cut(localism, breaks_qt$brks,include.lowest=TRUE))
ggplot(mergeddistricts) + geom_sf(aes(fill=localism_cat)) + scale_fill_brewer(palette = "YlOrRd",name="Localism") + theme(plot.margin=grid::unit(c(0,0,0,0), "mm"))
ggsave("localism.pdf",units="mm",width=150,height=150)


localism.Moran <-  moran.test(mergeddistricts$localism, listw)

sink("raw-localism-global-moran.txt")
cat(round(as.numeric(localism.Moran$estimate[1]),digits=3))
sink()





# Ok. We also want clustermaps for the raw data 

# culturalthreat 
  local.moran.cultthreat<-localmoran(as.vector(mergeddistricts$cultthreat),listw)
  #Copy classification (quadrant, divison by mean)
  mergeddistricts$local.moran.typ.cultthreat <- attr(local.moran.cultthreat,"quadr")$mean
  # and set to NA by insignifacnt LISA values 
  mergeddistricts$local.moran.typ.cultthreat[local.moran.cultthreat[,5] >=0.05] <- NA
  # Plot
  ggplot(mergeddistricts) + geom_sf(aes(fill=local.moran.typ.cultthreat)) + guides(fill=guide_legend(title="Cultural threat")) + theme(plot.margin=grid::unit(c(0,0,0,0), "mm"))
  ggsave("cultthreat-clustermap.pdf",units="mm",width=150,height=150)

# Islamophobia 
  local.moran.islamophobia<-localmoran(as.vector(mergeddistricts$islamophobia),listw)
  #Copy classification (quadrant, divison by mean)
  mergeddistricts$local.moran.typ.islamophobia <- attr(local.moran.islamophobia,"quadr")$mean
  # and set to NA by insignifacnt LISA values 
  mergeddistricts$local.moran.typ.islamophobia[local.moran.islamophobia[,5] >=0.05] <- NA
  # Plot
  ggplot(mergeddistricts) + geom_sf(aes(fill=local.moran.typ.islamophobia)) + guides(fill=guide_legend(title="Islamophobia")) + theme(plot.margin=grid::unit(c(0,0,0,0), "mm"))
  ggsave("islamophobia-clustermap.pdf",units="mm",width=150,height=150)

# Populism
  local.moran.populism<-localmoran(as.vector(mergeddistricts$populism),listw)
  #Copy classification (quadrant, divison by mean)
  mergeddistricts$local.moran.typ.populism <- attr(local.moran.populism,"quadr")$mean
  # and set to NA by insignifacnt LISA values 
  mergeddistricts$local.moran.typ.populism[local.moran.populism[,5] >=0.05] <- NA
  # Plot
  ggplot(mergeddistricts) + geom_sf(aes(fill=local.moran.typ.populism)) + guides(fill=guide_legend(title="Populism")) + theme(plot.margin=grid::unit(c(0,0,0,0), "mm"))
  ggsave("populism-clustermap.pdf",units="mm",width=150,height=150)


# Aggression/Submission
  local.moran.rwa<-localmoran(as.vector(mergeddistricts$rwa),listw)
  #Copy classification (quadrant, divison by mean)
  mergeddistricts$local.moran.typ.rwa <- attr(local.moran.rwa,"quadr")$mean
  # and set to NA by insignifacnt LISA values 
  mergeddistricts$local.moran.typ.rwa[local.moran.rwa[,5] >=0.05] <- NA
  # Plot
  ggplot(mergeddistricts) + geom_sf(aes(fill=local.moran.typ.rwa)) + guides(fill=guide_legend(title="Authoritarianism")) + theme(plot.margin=grid::unit(c(0,0,0,0), "mm"))
  ggsave("rwa-clustermap.pdf",units="mm",width=150,height=150)


# Regional exclusion 
  local.moran.regexclusion<-localmoran(as.vector(mergeddistricts$regexclusion),listw)
  #Copy classification (quadrant, divison by mean)
  mergeddistricts$local.moran.typ.regexclusion <- attr(local.moran.regexclusion,"quadr")$mean
  # and set to NA by insignifacnt LISA values 
  mergeddistricts$local.moran.typ.regexclusion[local.moran.regexclusion[,5] >=0.05] <- NA
  # Plot
  ggplot(mergeddistricts) + geom_sf(aes(fill=local.moran.typ.regexclusion)) + guides(fill=guide_legend(title="Place resentment")) + theme(plot.margin=grid::unit(c(0,0,0,0), "mm"))
  ggsave("regexclusion-clustermap.pdf",units="mm",width=150,height=150)


# Localism 
  local.moran.localism<-localmoran(as.vector(mergeddistricts$localism),listw)
  #Copy classification (quadrant, divison by mean)
  mergeddistricts$local.moran.typ.localism <- attr(local.moran.localism,"quadr")$mean
  # and set to NA by insignifacnt LISA values 
  mergeddistricts$local.moran.typ.localism[local.moran.localism[,5] >=0.05] <- NA
  # Plot
  ggplot(mergeddistricts) + geom_sf(aes(fill=local.moran.typ.localism)) + guides(fill=guide_legend(title="Localism")) + theme(plot.margin=grid::unit(c(0,0,0,0), "mm"))
  ggsave("localism-clustermap.pdf",units="mm",width=150,height=150)

# Now map/analyse the BLUPs (predicted random shocks)
# Join again to get BLUPs from Stata

blups <- haven::read_dta("district-blups.dta")
mergeddistricts <-inner_join(mergeddistricts,blups,by="GKZ")


# Start with cultural threat 
breaks_qt <- classIntervals(mergeddistricts$cultthreatblup, n = 7, style = "quantile")
mergeddistricts <- mutate(mergeddistricts, cultthreatblup_cat = cut(cultthreatblup, breaks_qt$brks))
ggplot(mergeddistricts) + geom_sf(aes(fill=cultthreatblup_cat)) + scale_fill_brewer(palette = "YlOrRd",name="Cultural threat") + theme(plot.margin=grid::unit(c(0,0,0,0), "mm"))
  ggsave("cultthreatblup.pdf",units="mm",width=150,height=150)

cultthreatblup.Moran <-  moran.test(mergeddistricts$cultthreatblup, listw)
sink("blup-cultthreatblup-global-moran.txt")
cat(round(as.numeric(cultthreatblup.Moran$estimate[1]),digits=3))
sink()

# Islamophobia 
breaks_qt <- classIntervals(mergeddistricts$islamophobiablup, n = 7, style = "quantile")
mergeddistricts <- mutate(mergeddistricts, islamophobiablup_cat = cut(islamophobiablup, breaks_qt$brks))
ggplot(mergeddistricts) + geom_sf(aes(fill=islamophobiablup_cat)) + scale_fill_brewer(palette = "YlOrRd",name="Islamophobia") + theme(plot.margin=grid::unit(c(0,0,0,0), "mm"))
  ggsave("islamophobiablup.pdf",units="mm",width=150,height=150)

islamophobiablup.Moran <-  moran.test(mergeddistricts$islamophobiablup, listw)
sink("blup-islamophobiablup-global-moran.txt")
cat(round(as.numeric(islamophobiablup.Moran$estimate[1]),digits=3))
sink()

# RWA
breaks_qt <- classIntervals(mergeddistricts$rwablup, n = 7, style = "quantile")
mergeddistricts <- mutate(mergeddistricts, rwablup_cat = cut(rwablup, breaks_qt$brks))
ggplot(mergeddistricts) + geom_sf(aes(fill=rwablup_cat)) + scale_fill_brewer(palette = "YlOrRd",name="Authoritarianism") + theme(plot.margin=grid::unit(c(0,0,0,0), "mm"))
  ggsave("rwablup.pdf",units="mm",width=150,height=150)

rwablup.Moran <-  moran.test(mergeddistricts$rwablup, listw)
sink("blup-rwablup-global-moran.txt")
cat(round(as.numeric(rwablup.Moran$estimate[1]),digits=3))
sink()


# Populism

breaks_qt <- classIntervals(mergeddistricts$populismblup, n = 7, style = "quantile")
mergeddistricts <- mutate(mergeddistricts, populismblup_cat = cut(populismblup, breaks_qt$brks))
ggplot(mergeddistricts) + geom_sf(aes(fill=populismblup_cat)) + scale_fill_brewer(palette = "YlOrRd",name="Populism") + theme(plot.margin=grid::unit(c(0,0,0,0), "mm"))
  ggsave("populismblup.pdf",units="mm",width=150,height=150)

populismblup.Moran <-  moran.test(mergeddistricts$populismblup, listw)
sink("blup-populismblup-global-moran.txt")
cat(round(as.numeric(populismblup.Moran$estimate[1]),digits=3))
sink()


# Regional exclusion / place resentment 
breaks_qt <- classIntervals(mergeddistricts$regexblup, n = 7, style = "quantile")
mergeddistricts <- mutate(mergeddistricts, regexblup_cat = cut(regexblup, breaks_qt$brks))
ggplot(mergeddistricts) + geom_sf(aes(fill=regexblup_cat)) + scale_fill_brewer(palette = "YlOrRd",name="Regional exclusion") + theme(plot.margin=grid::unit(c(0,0,0,0), "mm"))
  ggsave("regexblup.pdf",units="mm",width=150,height=150)

regexblup.Moran <-  moran.test(mergeddistricts$regexblup, listw)
sink("blup-regexblup-global-moran.txt")
cat(round(as.numeric(regexblup.Moran$estimate[1]),digits=3))
sink()

# Clustermaps 


# OK. Look at LISA, filter out insignifacnt values 
# Cultural threat 

local.moran.cultthreatblup<-localmoran(as.vector(mergeddistricts$cultthreatblup),listw)
#Copy classification (quadrant, divison by mean)
mergeddistricts$local.moran.typ.cultthreatblup <- attr(local.moran.cultthreatblup,"quadr")$mean
# and set to NA by insignifacnt LISA values 
mergeddistricts$local.moran.typ.cultthreatblup[local.moran.cultthreatblup[,5] >=0.05] <- NA
# Plot
ggplot(mergeddistricts) + geom_sf(aes(fill=local.moran.typ.cultthreatblup)) + guides(fill=guide_legend(title="Cultural threat")) + theme(plot.margin=grid::unit(c(0,0,0,0), "mm"))
ggsave("cultthreatblup-clustermap.pdf",units="mm",width=150,height=150)

# Islamophobia

local.moran.islamophobiablup<-localmoran(as.vector(mergeddistricts$islamophobiablup),listw)
#Copy classification (quadrant, divison by mean)
mergeddistricts$local.moran.typ.islamophobiablup <- attr(local.moran.islamophobiablup,"quadr")$mean
# and set to NA by insignifacnt LISA values 
mergeddistricts$local.moran.typ.islamophobiablup[local.moran.islamophobiablup[,5] >=0.05] <- NA
# Plot
ggplot(mergeddistricts) + geom_sf(aes(fill=local.moran.typ.islamophobiablup)) + guides(fill=guide_legend(title="Islamophobia")) + theme(plot.margin=grid::unit(c(0,0,0,0), "mm"))
ggsave("islamophobiablup-clustermap.pdf",units="mm",width=150,height=150)

# RWA 

local.moran.rwablup<-localmoran(as.vector(mergeddistricts$rwablup),listw)
#Copy classification (quadrant, divison by mean)
mergeddistricts$local.moran.typ.rwablup <- attr(local.moran.rwablup,"quadr")$mean
# and set to NA by insignifacnt LISA values 
mergeddistricts$local.moran.typ.rwablup[local.moran.rwablup[,5] >=0.05] <- NA
# Plot
ggplot(mergeddistricts) + geom_sf(aes(fill=local.moran.typ.rwablup)) + guides(fill=guide_legend(title="Authoritarianism")) + theme(plot.margin=grid::unit(c(0,0,0,0), "mm"))
ggsave("rwablup-clustermap.pdf",units="mm",width=150,height=150)

# Populism

local.moran.populismblup<-localmoran(as.vector(mergeddistricts$populismblup),listw)
#Copy classification (quadrant, divison by mean)
mergeddistricts$local.moran.typ.populismblup <- attr(local.moran.populismblup,"quadr")$mean
# and set to NA by insignifacnt LISA values 
mergeddistricts$local.moran.typ.populismblup[local.moran.populismblup[,5] >=0.05] <- NA
# Plot
ggplot(mergeddistricts) + geom_sf(aes(fill=local.moran.typ.populismblup)) + guides(fill=guide_legend(title="Populism")) + theme(plot.margin=grid::unit(c(0,0,0,0), "mm"))
ggsave("populismblup-clustermap.pdf",units="mm",width=150,height=150)


# Regional exclusion BLUP
local.moran.regexblup<-localmoran(as.vector(mergeddistricts$regexblup),listw)
#Copy classification (quadrant, divison by mean)
mergeddistricts$local.moran.typ.regexblup <- attr(local.moran.regexblup,"quadr")$mean
# and set to NA by insignifacnt LISA values 
mergeddistricts$local.moran.typ.regexblup[local.moran.regexblup[,5] >=0.05] <- NA
# Plot
ggplot(mergeddistricts) + geom_sf(aes(fill=local.moran.typ.regexblup)) + guides(fill=guide_legend(title="Regional exclusion")) + theme(plot.margin=grid::unit(c(0,0,0,0), "mm"))
ggsave("regexblup-clustermap.pdf",units="mm",width=150,height=150)

